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The emergence of vortex-bright soliton dipoles in two-component Bose-Einstein condensates 
through bifurcations from suitable eigenstates of the underlying linear system is examined. These 
dipoles can have their bright solitary structures be in phase (symmetric) or out of phase (anti- 
symmetric). The dynamical robustness of each of these two possibilities is considered and the 
out-of-phase case is found to exhibit an intriguing symmetry-breaking instability that can in turn 
lead to tunneling of the bright wavefunction between the two vortex "wells" . We interpret this phe- 
nomenon by virtue of a vortex-induced double well system, whose spontaneous symmetry breaking 
leads to asymmetric vortex-bright dipoles, in addition to the symmetric and anti-symmetric ones. 
The theoretical prediction of these states is corroborated by detailed numerical computations. 
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I. INTRODUCTION 

Within the booming field of Bose-Einstein condensation (BEC), the study of coherent nonlinear states has its own 
considerable history [IH1]. The original explorations were in the setting of repulsive interatomic interactions, and 
especially so in the context of one-component settings, starting over a decade ago [SHH]- These works chiefly focused 
on the dark soliton i.e., the prototypical nonlinear state therein, yet they were considerably hampered by instability 
effects induced by the dimensionality of the system and/or the presence of thermal effects. Nevertheless, more recent 
efforts using a variety of techniques have been far more successful in generating robust (dark) solitary wave states. 
Such techniques include phase- imprinting/density engineering |10H12j . matter- wave interference [131 114) . or dragging 
localized defects through the BECs [15] . The two-dimensional generalization of such (dark) states has as its prototype 
the vortex waveform, which became possible |16j by using a phase-imprinting method between two hyperfine spin 
states of a 87 Rb BEC [T7]. Subsequent efforts involved the stirring of the BECs [T5] above a certain critical angular 
speed [TM22] which, in turn, led to the production of few vortices [22] and even of very robust vortex lattices [23] . 
Other techniques including dragging obstacles through the BEC [H] or the nonlinear interference of condensate 
fragments 25 have been also used for the production of unit-charge vortices. Higher-charged vortex structures were 
produced j2Bj and their dynamical (in)stability has been examined. 

All of the above explorations were developed in the context of one-component BECs. Yet, solitary wave states also 
exist in multi-component BEC settings. In that context, of growing interest within the past few years has been the 
study of dark-bright (DB) solitons that are supported in two-component [27] and even spinor [28] condensates. These 
states can be thought of as "symbiotic" , in that the bright second component could not be sustained in the absence 
of the trapping dark first component. Robust such states were first observed in the experiment of Ref. [10] by means 
of a phase-imprinting method. 

This, in turn, has led to experimental studies of numerous features of these multi-component waves including the 
realization of DB soliton trains [53], DB soliton oscillations and interactions [30, 31 , as well as the possibility to create 
dark-dark breathing counterparts of these states (and multi-wave generalizations thereof) j32j|33]. In two dimensions, 
generalizations of these states have been proposed in the form of vortex-bright (VB) solitons, which were introduced 
about a decade ago (see e.g. [33] and references therein) and were recently further explored in the work of [55] . 

It is with the VB waveforms that we will concern ourselves in the present study. Early experiments, such as the 
one of [35], have illustrated the feasibility of realization of these states. Additionally, recent studies of dynamical 
phenomena in two-component condensates with considerable temporal and spatial resolution and control |37H39j 
suggest that the relevant coherent states can be explored further. Our aim here is to explore this potential beyond the 
level of a single vortex-bright soliton entity. In particular, recently in the one-component setting, we have examined 
the bifurcation of few- vortex clusters (see [30] and references therein), most notably the vortex dipole, but also the 
vortex tripole, quadrupole, vortex polygons and larger scale crystals |41) . It is then natural to expect that similar 
bifurcations will arise in the two-component setting and, in fact, that the bright solitons that are "trapped" within 
the vortex states will potentially bear different relative phases (e.g. in phase or out of phase, as is the case with 
DB solitons in [3T]). It is relevant to mention here that particular (in phase) realizations of such states have been 
very recently proposed as realizable by means of dragging numerical experiments in the immiscible regime of the 
pseudo-spinor system in |42j . The excitations resulting by dragging a laser beam through the two-component system 
had earlier been explored in the miscible regime in [43] . 

The vortex-bright dipoles, consisting of a vortex pair of opposite circulation and a corresponding trapped bright 
soliton pair, will be the main theme herein. We will start by providing the background of the relevant model and 
theoretical setup, as well as a brief review of the properties of a single VB soliton in section II. Then, in section III, we 
will explain the different types of [in phase (IP), or out of phase (OOP)] VB dipoles that exist and will present their 
bifurcation from the (nonlinear) continuation of underlying linear states in the form of dark-bright soliton stripes. In 
section IV, we examine the stability of the VB dipoles and recognize an instability of the OOP states and its symmetry 
breaking and tunneling implications. These are then theoretically explained in the form of an effective double well 
theory (as induced by the vortices on the bright components) and its symmetry breaking bifurcations, which gives 
rise to rather unexpected genuinely asymmetric VB dipole states. We summarize our findings in section V, where we 
also present some conclusions and possible directions for future studies. 

We should mention in passing that although our principal focus herein will stem from BEC and atomic physics 
considerations, relevant topics and ideas are, in principle, relevant for nonlinear optics as well. In particular, structures 
such as dark solitons [33] and vortices [35] [35] have been extensively studied in the latter field. In fact, DB-soliton 
states were also first observed in optics experiments, where they were created in photorefractive crystals [47], while 
their interactions were discussed in Ref. 48J . It is thus natural to expect that the combination of multi-component 
and multi-dimensional settings therein would also yield further potential for the realizability and observation of the 
coherent structures analyzed in the present work. 
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II. MODEL SETUP 



Our starting point will be the setting of quasi two-dimensional repulsive binary BECs, whose mean field description 
is given by the following set of equations: 



This coupled dimcnsionless set of Gross-Pitaevskii (GP) equations describes the time evolution of the two components' 
order parameters ipj, j £ {1,2}. Time, length, energy and densities \ipj\ 2 are measured in units of u~\ a z , hw z and 
(2\/27r\ai2\a z )~ 1 , respectively; oj z and a z denote the oscillator frequency and length in the frozen z-direction, while 
an, aii and a\i refer to the intra- and intercomponent scattering lengths. In the resulting dimensionless form of the 
equation, the coupling constants are gj — ajj/\a\2\, and <j\2 denotes the sign of a\2. 

In the following, we will analyze the relevant case of binary condensates composed of 87 Rb atoms in the two spin 
states (F = 1, mp = —1) and (F = 2, rap = 1), leading to dimensionless coupling constants of approximately 
g\ = 1.03, (72 = 0.97, <7i2 = +1 H], and we will exclusively consider isotropic harmonic potentials V(x,y) — Q 2 (x 2 + 
y 2 )/2 = Sl 2 r 2 /2, where 51 = w r /w z denotes the ratio of the in-plane and out-of-plane oscillator frequencies and is fixed 
to 0.2. This trap strength is only selected for reasons of computational convenience (and experimental realizability) , 
but the phenomenology presented below will not depend in any critical way on this selection for quasi-two-dimensional 
BECs. Stationary solutions are obtained by factorizing tpj(x,y,t) = exp(—ifijt)(f>j(x,y), j € {1,2}, where [L\ and /i2 
are the two components' chemical potentials. 

The prototypical stationary solution of Eq. ([!]) that will be the building block for our considerations is the single 
vortex-bright soliton state, whose density and phase profiles are shown in Fig. [T] For such a state, the first component 
supports a vortex, which acts as an effective potential well for the second component (of course, the alternative 
arrangement also exists where the role of the components is interchanged). However, we will focus solely on the 
former case due to its unconditional stability; as discussed e.g. in [341 135 ) . in the case where the components are 
interchanged, parametric regimes of instabilities may arise. 

The linear stability (so-called Bogolyubov-de Gennes or BdG) analysis is employed in order to consider the fate of 
small amplitude perturbations and the potential robustness of the solutions. This consists of imposing a perturbation 
to the stationary solutions <j)j above in the form: 



This leads at 0(s) (where e is a formal small parameter) to an eigenvalue problem for the (in principle, complex) 
frequency of excitations u> and the corresponding eigenvector [a(x,y),b(x,y),c(x,y),d(x,y)] T . Further details about 
the mathematical structure of the BdG two-component problem can be found in [34j . For our purposes, it suffices 
to note that the Hamiltonian structure of the resulting eigenvalue problem enforces that if w is an eigenfrequency so 
are uj* , —ui and — uj* . Hence, if eigenfrequencies with Im(w) =^ exist, then the solution is deemed to be dynamically 
unstable. 

A prototypical example of the BdG spectrum of a VB solitary wave is provided in the right panel of Fig. [T] In 
this example, the spectrum is offered as a function of the (rescaled) number of atoms in the bright soliton component 
N2 = J dxdy\ip2\ 2 , while the chemical potential of the dark (vortex) component stays fixed at [i\ — 5.2. One important 
observation to make here is that there is a single negative energy mode |34) in the spectrum of the vortex for small 
(or vanishing) N2 . This mode is illustrated by the square markers (which are red in the online version) in the relevant 
panel and is well-known to correspond to the precession of the vortex within the parabolic trap (see e.g. |40j and 
references therein). However, it is noteworthy that this mode decreases in frequency due to the presence of the bright 
component and ultimately crosses the zero frequency point. Yet, this crossing does not produce an instability; the 
relevant eigenfrequency pair remains real but now the energy of the mode is positive signaling the transition of the 
vortex state from a saddle in the energy landscape into a local minimum thereof (this is a setting analogous to what 
is observed for a vortex in the presence of rotation; see e.g. the relevant discussion of [IS]). It is important to note 
that this observation is in agreement with the results described in [36] , where it was observed that filled vortex cores 
exhibit slower precessions, and thus a decrease of the precession frequency is expected. It is relevant to note that 
similar observations have also been obtained in the case of dark-bright solitons, originally through the theoretical 
analysis of |27j and have been experimentally confirmed in the work of [30j . 



idtih(x,y,t) = --(dl+dl) + V{x,y)+g 1 \ip 1 \ 2 + <j 12 \ip 2 \ 2 i>x{x,y,t) 

id t Mx,y,t) = - l -{d 2 x + d 2 y ) + v{x lV ) + g2m 2 + o- 12 \^\ 2 ihfav,t). 



(1) 




(2) 



(3) 
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FIG. 1. Left panel: numerically calculated vortex-bright soliton solution to the two-component Gross-Pitaevskii equation 
obtained for Ni — 2000 and N2 = 100. The density (left) and the phase (middle) of the two components is shown. Right 
panel: BdG spectrum of the vortex-bright soliton for hi = 5.2 as a function of iVjj. Just the real part of the eigenfrequencies 
u> is shown, as the imaginary part is found to be identically zero for all values of N2. Negative energy modes are indicated by 
darker square markers (red in the online version). 



III. VORTEX-BRIGHT SOLITON DIPOLE STATES 



Using the above fundamental building blocks, namely the single VB solitons, we now look for bound states containing 
multiple such entities. In the same spirit as in the work of [40] (see also e.g. [50] and references therein), the prototypical 
relevant bound state is the VB dipole. We have been able to identify two such states, which are both shown in Fig. [2] 
for Ni = 70 and N2 = 20. In both cases the first component contains two vortices located symmetrically with respect 
to the trap center, which are filled by the second component. The difference between the two dipoles is evident in the 
phase profile: the two bright solitons filling the vortex cores can have either the same phase, or there may be a phase 
difference of n between them. 

The emergence of the dipole state branches can be thought of in terms of a bifurcation picture in an equivalent 
way to the one described for one-component BECs in [3D]. The two possible dipole branches, distinguished from the 
presence or absence of the phase jump of 7r in the bright component, arise at critical values of N\ from two states, 
which can be regarded as generalizations of the dark soliton stripe observed in one-component condensates [4]. The 
density and phase profiles of these two solitonic states are shown in Fig. [3] where it is evident that the two states are 
again distinguished by the presence or absence of a 7r-phase jump in the bright component. 

In the first case of the IP VB dipole, the parent state can be seen to consist of a dark soliton stripe in the first compo- 
nent accompanied by a bright stripe in the second component, which is trapped inside the dark stripe. In the notation 
of the corresponding linear limit of the two single-particle Hamiltonians, the relevant state of the first component is the 
|1, 0) state (i.e., the first excited state in the x-direction) while that of the second component is the |0, 0) ground state 
of the latter. More generally, the symbolism \m, n) is used to denote the states H m (\^x)H n (\/Tly)e^ nr / 2 where H m 
stands for the m-th Hermite polynomial and normalization constants have been omitted. These are stationary states 
of the underlying linear problem with eigenvalue /1 = £l(m + n + 1), from which the soliton stripe solutions arise in the 
presence of the effective nonlinearity induced by interatomic interactions. It is the nonlinearity and immiscibility of 
the two components that leads to the |0, 0) state of the second component being elongated (as opposed to circularly 
symmetric) in the dark-bright stripe of the top left panel of Fig. [3] Corresponding states in a ring (as well as in 
diagonal stripes) form have recently been addressed in [5TJ. In the case of the top right panel of Fig. [3] and the OOP 
VB dipole, the parent state is the |1,0) of the first component coupled to the |0, 1) of the second component. 



5 




FIG. 2. Left panels: Density (first column) and phase (second column) of an in-phase vortex-bright dipole. Right panels: 
Density (third column) and phase (fourth column) of the out-of-phase vortex-bright dipole. Both examples are for Ni = 70 
and N 2 = 20. 



While the bright component remains roughly invariant as the transition point is crossed, the dark component devel- 
ops two singular points in the phase corresponding to the density vanishings and the opposite circulation characteristic 
of a vortex dipole 40 . The relevant bifurcations occur approximately at the same critical N\ (iVf ss 40, with 
fixed at 20). This is explicitly depicted in the middle panels of Fig. [3] for the two different "stripe" states, where the 
bifurcation diagrams which determine the relation Ni(ni) are shown. These approximately equal critical values of iVi 
can, arguably, be expected due to the small number of particles of the second component, a statement which becomes 
exact in the limit of N2 — > 0. 

The BdG analysis of the two stripe states is shown in the bottom panel of Fig. [3] and allows to relate the emergence 
of the dipole state to the changing of stability of the parental branch (the IP and OOP solitonic states). 

From the BdG spectra, in fact, one observes the emergence of an imaginary eigenfrequency pair for Nf r ~ 40, which 
is the value of Ni at which the bifurcation actually occurs. Imaginary modes are here depicted by circle (pink in the 
online version) markers. The "*" (green in the online version) markers, are used to identify complex modes which arise 
from the collision of two eigenmodes with opposite Krein sign. Such weak instabilities emerging in correspondence 
with collisions of two real modes are often observed in the numerically obtained spectra in this and several other 
works 40J. However, they are usually not considered, as the parametric regimes where thay appear are very small 
and the observed instabilities are weaker than the ones associated with purely imaginary eigenfrequencies. 

Let us also note in passing that the IP stripe bears only one negative energy mode (which can be seen to lead 
to oscillatory instabilities), while the OOP one bears two such anomalous modes (which, as for the single vortex 
bright soliton, are indicated by square, red in the online version, markers). This is intuitively expected by the out- 
of-phase or excited state nature that the former bears only in the first component, while the latter has that type 
of structure in both components. I.e., in the second component, the former configuration features a ground state, 
symmetric waveform, while the latter has an anti-symmetric first excited state waveform. This type of characteristics 
will be of critical relevance to the considerations that follow below, regarding tunneling effects and symmetry-breaking 
bifurcations. 



G 





FIG. 3. Top panels: density and phase of the in-phase dark-bright soliton stripe (first and second column) and of the out- 
of-phase corresponding state (third and fourth column), both at Ni — 38, N2 = 20. The middle panels show, respectively, 
the bifurcation of the vortex-bright soliton dipole branch from the corresponding parent states (of the top panels). Finally, 
the bottom panel contains the BdG analysis of the parent states (on the left and right sides we see the IP and the OOP 
case, respectively) of the top panels, illustrating the instability (through an imaginary eigenfrequency, signaled by pink circle 
markers) that leads to the emergence of the vortex-bright dipoles. In addition, for the out-of-phase case, the lines formed by 
"* " (green in the online version) markers correspond to the existence of oscillatory instabilities through complex eigenfrequency 
quartets. In this case there are two anomalous modes (depicted by red square markers), while there is only one such for the in 
phase stripe. 



7 



IV. TUNNELING DYNAMICS AND SYMMETRY-BREAKING BIFURCATIONS 



Let us now come back to the vortex-bright soliton dipoles and especially their dynamical stability through the BdG 
analysis. The in-phase VB dipole appears to be stable for arbitrary values of iVi and N2, while the out-of-phase VB 
dipole shows purely imaginary modes arising at certain values of Ni and 7V 2 (when the state is scanned over Ni or 
N2, respectively). Exciting this unstable mode for the OOP dipole with Ni = 220 and N2 — 1 by adding some white 
noise and letting the perturbed states propagate in time, the evolution shown in Fig. [4] is obtained. 




FIG. 4. Time propagation of the perturbed out-of-phase vortex-bright soliton dipole state at Ni = 220 and A^2 = 1. The top 
and bottom lines depict the density profiles of the dark and bright component, respectively, evaluated at five moments of the 
simulated evolution. The elapsed (dimensionless) time is indicated in the panels' upper left corners. The substantial addition 
of white noise (that seeds the instability) can be seen in the pixelization of the corresponding dynamical panels. 



The time evolution of the density profiles shows the size of one of the bright solitons increasing while the other one 
is getting smaller, then the opposite situation is observed and this trend with enlarging and shrinking bright solitons 
is repeating in the course of the propagation. After an initial transient phase during which the bright component's 
asymmetry builds up (first to second panel of Fig. [4]), we find periodic oscillations. This observed time evolution has 
an immediate physical interpretation as concerns the bright component. A fraction of the corresponding particles 
within the second component move from one vortex core to the other during the time propagation of the system. 
This can be explained by considering the dark component's density as acting in the form of an effective potential 
for particles of the bright one. I.e., the particles of the second component can be considered as tunneling within the 
vortex-induced (i.e., formed by the vortex cores) double well potential. Within this very potential, this suggests the 
possibility of a spontaneous symmetry breaking bifurcation, as responsible for the observed phenomenology. 

In order to validate this assumption, we consider the analytical model, developed in |52j . based on a two-mode 
expansion, which is used to determine stationary states and to study the symmetry-breaking bifurcations occuring 
in one-dimensional, single-component repulsive BECs confined in double-well potentials. The model predicts that for 
a two-mode decomposition of the equation involving a symmetric ground state and an anti-symmetric first excited 
state, there exists a symmetry breaking bifurcation both for attractive and for repulsive interactions. The bifurcation 
emerges from (and destabilizes) the anti-symmetric, first-excited state in our repulsive interaction case (while it stems 
from and destabilizes the symmetric ground state in the attractive interaction case). This is consonant with the 
principal observation above of the generic stability of the IP VB, which has a symmetric second component and the 
symmetry breaking associated destabilization of the OOP VB with the anti-symmetric second component. 

The key quantitative observation now is that if we freeze the first component (assuming that it forms the vortex-core 
double well potential), then the theory of [52] can be applied directly but for the effective double well potential in the 
form 

V e s(x,y) — V{x,y) + \ipi{x,y)\ 2 , where V(x,y) is the trapping harmonic potential. The equation for the second 
component can then be rewritten as 



id t tp 2 (x,y,t) = 



~Ia 2D + V eS (x,y) + salVaf 



tp 2 (x,y,t) 



(4) 
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which upon the rescaling ip2(x,y) = y/g2ip2(x,y) leads to a standard single-component (two-dimensional) double well 
setting. It should be noted that the potential of the double well clearly depends on Ni (through its dependence on 
\i^i(x, y)\ 2 ). For each N% (in multiples of 5) in the interval [0, 300] the density profile of the dark component of the VB 
dipole state with N2 = has been utilized (to form the corresponding effective potential). For each considered value 
of N±, one can diagonalize the obtained Hamiltonian H = — \£±2D + Vcfi, and keep the two lowest energy eigenstates, 
namely the symmetric ground state and the anti-symmetric first excited state of the single particle operator with the 
effective potential. These, denoted hereafter as u (x,y) and Ui(x,y), will be used for the two mode reduction in the 
spirit of [52] , 

The fundamental analytical prediction of the two-mode approximation is that as N2 (and [12) is increased, a critical 
point will be reached, given by 



N ** = (W^A l)92 > (5) 

whereafter the anti-symmetric branch (interpreted here as the OOP VB dipole) will become unstable. Past this 
critical point, an asymmetric branch (i.e., an asymmetric VB dipole) should emerge as a stable configuration. In the 
expression of Eq. (15]) A\ = J ufdxdy and B = f UQU 2 dxdy are overlap integrals of the two lowest energy eigenstates 
of the linear Schrodinger problem and Alu — uji — luq, where w and ui\ are the energy eigenvalues corresponding to u 
and u\. The presence of the factor (72 in the denominator ensures consistency with the scaling of %p2 discussed above. 

The critical chemical potential is calculated from the particle number by making use of the expression [i2cr = 
u + 3BN 2cr [52] . 

Repeating this procedure for each considered Ni, we get the critical values of /Z2, which are compared with the 
numerically obtained data; see the top right panel of Fig. [5j These latter are obtained by making a scan within 
a suitable interval of N2 of the OOP VB dipole, and then performing the BdG analysis of the resulting branch of 
solutions. In the ensuing BdG spectrum, the bifurcation point can be identified as the onset of instability of the 
anti-symmetric VB dipole. For a relevant example, see the top left panel of Fig. [5] for the scan obtained for N± = 220, 
where we can see the imaginary mode emerging at N2 ~ 0.6 (generally, the relevant critical values iVf are of order 
unity). At this critical value of N2, for which the antisymmetric dipole becomes unstable (or, equivalently, at the 
corresponding critical value of ^2), we expect the bifurcation of a stable asymmetric state, which has been actually 
verified, as will be shown below. 

From the top right panel of Fig. [5j it is evident that the agreement between analytical predictions and numerical 
results is very good for what concerns fi2cr for Ni > 120. For lower values of N\ numerics and analytical predictions 
start to be in disagreement with each other and this has the following explanation: by decreasing N%, the height 
of barrier between the two wells forming the effective potential V e f[(x,y) — V(x,y) + \ipi(x, y)\ 2 , is also decreased. 
Thus, for low N%, the density \ijji(x,y)\ 2 provides a much lower contribution to the potential, and the effect of the 
harmonic trap is much more appreciable. This influences the form of the energy spectrum: the gap separating the 
almost equal lowest energy eigenvalues ojq and wi from the larger eigenvalues UJ2, w 3---, which is typically large for 
a high barrier, if the latter is lowered becomes smaller and smaller, until the spectrum of the harmonic oscillator, 
formed by equidistant eigenvalues, is reached. But, if this is the case, the two- mode approximation that one makes 
by projecting the problem on the eigenmodes uq and u\ fails as the contribution of the other eigenstates cannot be 
neglected anymore. 

Furthermore, another intuitive argument can be provided to explain the disagreement between predicted and 
numerical results observed with lowering N\: as was already said, the effective potential V c g(x, y) = V(x, y)+\ipi(x, y)\ 2 
has been calculated deriving the term \tpi{x, y)\ 2 from the density profile of the dark component of the VB dipole 
state with N 2 = and is considered to be fixed in the whole calculation. Our model looks for the bifurcation to occur 
at a finite value of iVjj, but neglects that, for this value, in the actual physical system, the density profile for the dark 
component is different from the one with N2 = 0. Thus, in doing this approximation, we do not take into account the 
effect of the second component on the first due to the interaction between them. Considering V e s independent of N2 
is expected to be valid if N% is large enough, as the dark component will then have a robust configuration with respect 
to variations due to the intercomponent interaction. But, decreasing Ni, varying N2 will start to influence both 
components in a sensitive way and, in particular, an evident variation of the density profile of the dark component is 
expected. Therefore, for low N±, the assumption that |V'i(s,2/)| 2 is independent of N2 is no longer appropriate and 
the model described above cannot be expected to be valid anymore. Therefore, in the top right panel of Fig. [3J values 
of Ni smaller than 100 are not taken into account. 

Let us now come back to the expected bifurcation from the anti-symmetric branch, once its destabilization occurs, 
towards an asymmetric VB dipole solution. Indeed, an example of this bifurcation has been illustrated in the bottom 
left panel of Fig. [5j The bottom right panel of the figure illustrates a prototypical example of the daughter state in 
the form of an asymmetric vortex-bright dipole. 



9 



0.25 




-10 10 -10 10 



x x 



FIG. 5. The BdG analysis of an antisymmetric VB out of phase dipole for Ni — 220 is shown in the top panel, as a function 
of N2. According to our convention, pink circle markers are related to unstable imaginary modes, while red square markers 
identify anomalous modes. The critical points for this instability (p2cr for different values of Ni) are given in the top right 
panel (red circle points) and are compared to the corresponding values of the theoretical two-mode prediction (solid black line). 
The bottom left panel shows the bifurcation diagram that demonstrates the existence of an asymmetric vortex-bright dipole, 
as bifurcating from the out-of-phase one, and the bottom right panel shows an example of the density and phase of such a 
symmetry-broken state. 
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The robustness of such an asymmetric VB dipole, combined with the Hamiltonian dynamics of the system, is 
what gives rise to the tunneling observations of Fig. [4] (since the anti-symmetric solution becomes a saddle, and 
the asymmetric one a center around which the dynamics revolves). Naturally, there are two such asymmetric states, 
depending on which vortex-induced well picks up part of the bright mass of the other, justifying the pitchfork character 
of the bifurcation. The BdG spectrum of the asymmetric state (not shown here) reveals the complete stability of the 
latter, in accordance with the expectations from the above bifurcation theoretic arguments. 

So far, we have concentrated the discussion of the tunneling dynamics on the regime of small N2 , close the linear 
limit of the effective Hamiltonian introduced by assuming the dark component to be frozen. While this is essential 
for the simplified semi-analytical model introduced above to make the connection to symmetry-breaking in a double 
well potential, we have found direct numerical evidence that the out-of-phase VB dipole is unstable towards tunneling 
dynamics in the bright component even in parameter regimes where N2 is considerably larger, see Fig. [6] for an 
example. In this case the larger population of the bright component now leads to a stronger back-coupling to the 
dark component, and it can be observed that the vortex cores forming the dipole in the dark component change size 
following the trend of the bright solitons filling them. 

While in Fig. [4j the tunneling oscillations, past the initial transient of asymmetry build-up, directly turn periodic, 
this is no longer the case here. Instead, in the run shown in Fig. [6j the bright component first almost completely 
tunnels to the lower well, but then oscillates between this asymmetric and a more symmetric occupation of the two 
wells (see the first four timesteps shown). Only after that, a majority of bright component atoms enters the upper 
well (second to last timestep), and then directly oscillates back to approximately equal occupation of the vortex cores 
(last timestep). The conclusion of this first period is followed by a series of nearly equal period similar oscillation 
steps. 
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FIG. 6. Time propagation of an out-of-phase vortex-bright soliton dipole, perturbed with white noise at Ni = 390 and N2 = 50. 

Another relevant difference to Fig. [4] is that the tunneling dynamics is accompanied by a rotational motion of the 
vortices, while for smaller 7V 2 the dipole stayed aligned along the trap's axis during the propagation. This may be 
taken as a hint that for a larger JV2/JV1 ratio the interspecies interaction may induce a coupling between the unstable 
tunneling-type mode of the bright component and the precessional degrees of freedom of the vortices in the dark 
component. 



V. CONCLUSIONS AND FUTURE CHALLENGES 

In the present work, we have revisited the vortex-bright solitary wave states and have illustrated their prototypical 
generalization to vortex-bright soliton dipoles. This is a first step towards the realization of two-component clusters 
of such entities. In the process, we have explored a number of interesting features. Firstly, even a single vortex-bright 
soliton was found to exhibit a transition from a saddle point in the energy to a local minimum, as its bright component 
becomes more significant. Secondly, two types of VB dipoles were identified, in analogy with their dark-bright one- 
dimensional counterparts. The first of them had the bright pulses (trapped inside the vortex cores) be in phase, 
while the other one had them as out of phase. We have identified the bifurcation of these states, as stemming from 
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the corresponding dark-bright stripe states, one of which has the bright second component in phase (stemming from 
the ground symmetric state) and one of which has it out of phase (stemming from the first excited anti-symmetric 
state). Subsequently, the dynamical stability of the daughter VB dipole states emerging from these bifurcations 
was quantified. It was found that the in-phase one is generically stable, while the out of phase structure suffers an 
exponential instability (past a relevant critical point). The latter gave rise to tunneling dynamics and spontaneous 
symmetry breaking manifestations. This phenomenology was understood on the basis of an effective double well 
model, where the vortex cores played the role of the wells and led to asymmetric VB dipoles. 

This study paves the way for a more detailed understanding of multi-component structures. One of the key items 
that remain open concerns the effective description of such states. In particular, it is remarkable that these states have 
both a solitonic character through their bright component and a vortex character through their dark one. It is then 
particularly interesting to examine how the effective equations characterizing the interactions of such entities look 
and what the corresponding interplay is between the vortex and the solitonic character. Another natural direction 
concerns generalizing the ideas presented herein to a setting with more vortex-bright dipoles e.g. three such, which 
might naturally be observable in two-component generalizations of the experiments of [53] . In fact, careful inspection 
of the figures presented herein (see e.g. the second pair of imaginary cigenfrequencies stemming from the stripe 
states' BdG analysis in the bottom panel of Fig. [3]) already suggests the bifurcation of such states. Furthermore, 
extending the present considerations to three-dimensions and towards a more detailed understanding of bound states 
of vortex rings |54j (and two-component generalization thereof) would also constitute an important theme for future 
explorations. Such studies are presently in progress and will be reported in future works. 
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